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Abstract 

Monte Carlo simulations of the 4-dimensional compact U(l) lattice gauge 
theory in the neighborhood of the transition point are made difficult by the 
suppression of tunneling between the phases, which becomes very strong as 
soon as the volume of the lattice grows to any appreciable size. This problem 
can be avoided by making the monopole coupling a dynamical variable. In 
this manner one can circumvent the tunneling barrier by effectively riding on 
top of the peaks in the energy distribution which meet for sufficiently large 
monopole coupling. Here we present an efficient method for determining the 
parameters needed for this procedure, which can thus be implemented at low 
computational cost also on large lattices. This is particularly important for a 
reliable determination of the transition point. We demonstrate the working of 
our method on a 16 4 lattice. We obtain an equidistribution of configurations 
across the phase transition even for such a relatively large lattice size. 



1. Introduction 



Recent investigations [1], |2|] of the compact U(l) lattice gauge theory in 4 dimen- 
sions have produced energy histograms indicative of a first-order transition. In the 
corresponding Monte Carlo simulations, however, the tunneling between the phases 
is strongly suppressed. In order to overcome the difficulties due to the lack of tran- 
sitions the authors of (JT] introduce an iterative reweighting for different (3, while the 
authors of || use a matching of hot and cold start results. The problem is that 
on larger lattices conventional algorithms are not able to induce transitions at all. 
Therefore, we have looked for algorithmic alternatives. 

To reduce the slowing down of the Monte Carlo algorithm in systems with a rough 
free-energy landscape the method of simulated tempering Q has been proposed and 
applied to the random-field Ising model. In this method the inverse temperature (3 is 
promoted to the status of dynamical variable, taking values which range over some 
definite set. In this manner one tries to utilize the fact that at lower (3 the free-energy 
barriers are lower. In an application of this procedure to spin-glass simulations f|] 
it has turned out, however, that adjusting the set of temperatures and handling the 
corresponding probability distribution in an efficient way is far from straightforward. 
Nevertheless it has been possible to develop a procedure |4j] leading to a reduction of 
slowing down comparable to the one obtained with the multicanonical method 0, 
and with the additional advantages of allowing full vectorization of the code and of 
providing the canonical ensemble directly. 

In the Potts model in 2 dimensions, the strength of the first-order transition de- 
creases with the number of the degrees of freedom q of the spins, the transition 
becoming of second order for q < 5. This has been used to set up an algorithm 
H in which q becomes a dynamical variable: by opening the easier pathway along 
the mountains of the joint probability distribution of q and energy, one avoids the 
need of relying, for large q, on the strongly suppressed tunneling for equilibrating 
the configurations. To implement transitions between different q cluster steps [|7| 
have been inserted. It turns out that by this algorithm one gains large factors in 
the autocorrelation times also in comparison to the multicanonical algorithm ||. 

Proceeding along similar lines, we have obtained an efficient algorithm for the U(l) 
gauge theory ||. We start from the Wilson action supplemented by a monopole term 
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S = (3 Y, (l-cose^J+A^lM^I, (1.1) 

fl>V,X p,x 



where M px is the monopole content of 3 dimensional cubes 11 . One finds that 



the strength of the first order transition decreases with A, the transition ultimately 
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becoming of second order. Thus, by making A a dynamical variable, we can again 
dispose of the tunneling transitions and proceed instead along the much easier path- 
way running over the top of the joint probability distribution P(E, A), E being the 
average plaquette energy. With the use of appropriate Metropolis steps for the link 
variables as well as A, moreover, one can make the algorithm fully vectorizable and 
parallelizable. 

Before running the dynamical-parameter algorithm one has to determine the posi- 
tion of the phase transition as function of A and some parameters in the generalized 
action, which serve to enforce the prescribed A distribution. On lattices of moderate 
size (e.g. 8 4 ) this is relatively easy because there is still some overlap between the 
peaks. On larger lattices determining these quantities is much more difficult and 
it becomes then crucial to perform the calculation without excessive computational 
cost. In the present paper we develop a method to achieve this goal. We demon- 
strate its effectiveness illustrating results obtained for a 16 4 lattice. We will see that 
our method enables us to observe transitions also on large lattices, which is very 
important for a reliable determination of the transition point. 

In Section 2 we outline the general features of the dynamical-parameter method. 
In Section 3 we derive relations among transition probabilities on which we will base 
the determination of the quantities required for the implementation of the algorithm. 
The detailed procedure followed for their calculation is described in Section 4. In 
Section 5 we will present some numerical results. 

2. Outline of method 

Conventional methods simulate the probability distribution 



where A is a fixed parameter. In order to make A a dynamical variable we consider 
H\{&) as the conditioned probability to get a configuration O given a value of A and 
prescribe a probability distribution /(A) to get the joint probability distribution 
//(©, A) = f(X)fj,\(Q). To simulate //(©, A) we need it in the form 




(2.1) 



/i(0,A) =exp(-S(e,\))/Z 



(2.2) 



where 



S(Q,\) = S x (Q)+g(X) 



(2.3) 
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and 

g(X) = -\og(f(X)Z/Z x ). (2.4) 

Eventually we will require that the values of A be approximately equidistributed, 
i.e. /(A) ~ const, which then gives g(X) ~ \ogZ\+ const. 

In our application of the algorithm each update of the link variables G MjX is followed 
by an update of A, which can take values from a discrete, ordered set X q with 
q = 1, . . . ,n. The individual update steps are Metropolis steps in both cases. For 
the A update we use a proposal matrix ^(8 q+ i iq i + 5 qs >+i + <5 g ,i<V,i + ^,nV,«) an d 
an acceptance probability min(l, exp(S'(9, X q ) — S{Q, X q i))). The above form of the 
proposal matrix implies that, if the current value X q is not extremal, then we choose 
as new candidate value for A one of the two neighboring values, A 9 _i or A 9+ i, with 
equal probability, whereas, if X q lies at the boundary of the set of possible values, we 
preselect either its (only) neighboring value or X q itself, again with equal probability. 

In order to implement the simulation, one must fix j3(X q ) (cfr. ( |1.1|) ) and g(X q ) 
(cfr. (|2.2| ) and (|2.3|)), for all values of q. We demand (3(X q ) ~ j3 w (X q ), where (3 W is 
the value of (3 which makes both phases equally probable. Our condition for fixing 
g(X q ) is /(A) ~ const. In order to determine (3(X q ) and g(X q ), we will use the fact 
that in a simulation the transition probabilities between neighboring values of A are 
very sensitive to these quantities. 



3. Transition probabilities 

To derive relations which can be used for the envisaged determination of /3(X q ) 
and g(X q ) we use the probability for the transition from a value X q to a neighboring 
value X q i 

W(Q, q; q') = ± min(l, exp(5(6, X q ) - S(Q, A,,)) (3. 1) 

and note that detailed balance implies 

/(A ff _0/i Ae _ 1 (e)^(e, q-i;q) = f(x q )» Xq (e)w{e, q - q - 1) . (3.2) 

Let us consider subsets K(q) of configurations with probability distributions 
proportional to fi\ q (0) and weight Wxiq) = 12eeK(q) A*a<,(0)- If we introduce the 
average transition probability for the set K(q) 

p K (q; q') = E ^,(6)1^(6, q; q') , (3.3) 
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by averaging ( |3.2| ) we obtain 



w K (q - 1) p K (q -l;q) = f(X g ) w K (q) p K (q] q ~ 1) ■ 



(3.4) 



We now apply ( |3.4j) to sets K c and Kh of configurations in the cold phase and 
in the hot phase separately. Because we are interested in cases where transitions 
between the phases are extremely rare, in practice it is easy to obtain sets of this 
type with numbers of configurations sufficient for the present purpose. Also, for the 
same reason, the corresponding equations can be considered to be independent. We 
assume that our two conditions, /3(A) = f3 w {X) and /(A) =const, are satisfied. The 
condition on f3 implies that the two phases are equally populated, i.e. wk c = WK h - 
Moreover, since the two subsets and K c essentially exhaust the whole set of 
configurations (in the cases we are considering the overlap is extremely small), all 
of the weights are, to a very good approximation, equal to |. Using this fact, the 
constancy of /(A), and ( ft.4| ) for K = K c and K = K h separately, we get a pair of 
equations which simplifies to 



This is what we will exploit to determine /3(X q ) and g(X q ). 

Our strategy will be to adjust P(X q ) and g(X q ), for known /?(A g _i) and g(\ q _i), in 
such a way that ( |3.5| ) holds. Starting from given /3(Xx) and arbitrarily chosen g(Ai), 
in this manner we can obtain (3{X q ) and g{X q ) for q — 2, . . . , n. 



To begin our procedure we select a value for Ai in the region where the peaks of 
the probability distribution associated to the two phases strongly overlap so that 
tunneling occurs frequently and /3(Ai) can easily be obtained by a conventional 
simulation. Because only the differences g(X q -i) — g{X q ) are relevant we can choose 
g(X\) arbitrarily. Then for q = 2, . . . , n we consecutively determine (3{X q ) and g(X q ) 
for known /3(A ? _i) and g(X q -i). 

In order to proceed from q — 1 to q we choose a new X q , approximately at the 
same distance from A 9 _i as in the previous steps. As a first rough approximation 
we obtain /3(A ? ) by extrapolation from the former values. At this point we use the 
sets of cold and hot configurations K c (q — 1) and K h (q — 1) at A ? _i, available from 



PK c (q-l;q) =PK B (q;q 

PK h (q-i;q) =PK h (q;q 



i) 
i) 



(3.5) 



4. Determination of /5(A) and g(X) 
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the previous step, and generate two new sets of 9 configurations K c (q) and Kh{q) 
at X q by short Monte Carlo runs with cold and hot start, respectively. For each set 
Ki(q') we can easily calculate the quantity 

PK t tf;<f) = -j^— E WW;?"), (4.i) 

where Nx^q 1 ) is the number of configurations in the set and W(Q, q'; q") is given by 
(|3.1| ) (of course, the variables q', q" stand for q — 1, q or q, q — 1, as appropriate). 

Since p approximates Q3.3|), this allows us to calculate the quantities px t which, for 
the correct choice of /3(X q ) and g(X q ), should satisfy (|3~5| ). We adjust then f3(X q ) and 
g(\q) until the equations (|3.5| ) are satisfied. In practice this takes only a very small 
amount of computer time. We obtain good estimates for (5{X q ) and g(X q ) though 
only approximate quantities enter ( |3.5|) because the peaks related to the phases 
vary only little with (3. In addition, the quantities pKi{q'] q") are used to adjust the 
distances between neighboring A values in such a way that one has roughly equal 
transition probabilities for all steps. 

After a larger number of q steps the errors may accumulate. Therefore we perform 
short runs of the dynamical-parameter algorithm to test whether the simulation 
does indeed travel along the mountains of the distribution in the hot as well as in 
the cold phase. If it gets stuck we slightly increase or decrease the couplings (3(X q ) 
in the region of A where the transitions fail. We determine then the corresponding 
values g(X q ) from the conditions 

Pk c (q- i;?) +PK h (q - i;g) =pk„{q',q- 1 )+PK h (q;q- x ) ( 4 - 2 ) 

and run the dynamical algorithm again. Typically one or two trials are sufficient. 

After performing the simulations with dynamical A, improved /3(\ q ) can be ob- 
tained by reweighting |T2[ the distribution at the values of A where deviations from 



the equidistribution of configurations in the cold and hot phase are seen to occur. 
Corresponding new values for g(X q ) are then obtained from ( j4.2| ). Alternatively 
improved values for g(X q ) can be obtained by replacing the current values with 
g(X q ) + \n(f(X q )). 



5. Numerical results 

Our method has made it possible to determine the phase transition region for a 
lattice as large as 16 4 using only a moderate amount of computer time. We have used 
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approximately 2 x 10 4 sweeps per each value of A to get the sets of configurations K c 
and Kh and approximately 4 x 10 4 sweeps, in total, in the short test runs with the 
dynamical algorithm. These preliminary calculations have been used to determine 
X q , (3{X q ) and g(X q ) following the procedure described in Section 4. Our results for 
these parameters are reproduced in Table 1. Altogether we used 25 values for A 
ranging from A = 0.4 (our starting point) down to A = 0. In our simulations with 
dynamical A we performed approximately 10 6 sweeps of the lattice and we observed 
a large number of transitions between the phases also on the 16 4 lattice. 

We define as location (3c of the phase transition the maximum of the specific heat. 
We have used reweighting [I2| in order to explore a range of (3 in the neighborhood 
of the value /3(X q ). As a matter of fact reweighting is necessary not only to find /3c, 
but also to determine accurately (3 W (the value of (3 where the configurations are 
equidistributed between the phases) since in order to implement the procedure of 
Section 4 we only needed to make the areas under the peaks approximately equal. 

Figure 1 shows (3c as function of A for the 16 4 lattice and also our earlier results 
H for the 8 4 lattice. In particular, for the 16 4 lattice at A = 0, we obtain the 
value (3c = 1.01084(5) , where the error has been estimated from the fluctuation of 
different samples. This confirms the result (3c = 1.01082(6) obtained in Ref. || by 
a matching procedure. 

In Figure 2 we show the distribution P(E, A) at (3 W (for A = 0.6 at (3c) which we 
got (after reweighting) from our simulations. The data have been obtained with the 
dynamical-parameter algorithm, except for A = 0.5 and A = 0.6 where the peaks 
overlap and the conventional Metropolis algorithm is adequate. Comparing with 
the corresponding figure for the 8 4 lattice in || the much stronger suppression of 
tunneling in the transition region is obvious. In fact, on the 8 4 lattice, because of 
the overlap between the peaks, there is still substantial tunneling. (For this reason, 
in our earlier simulations on the 8 4 lattice we could determine g(X q ) following the 
less sophisticated procedure based on Q4.2Q .) 

In regards to the efficiency of our algorithm versus conventional methods, the 
number of sweeps required to observe comparable numbers of tunnelings is greatly 
reduced already on an 8 4 lattice. One must make here a distinction (cfr. the 
discussion in 0) according to whether one is interested in all values of A, in which 
case our method produces all of the results in one stroke, or in a single A. In the 
latter case, since our method requires that one still simulates a whole range of A 
values, fairness requires that the the observed mean time between tunnelings be 
multiplied roughly by the number of A values considered. Even in this case, with an 
8 4 lattice there is still considerable gain, for example, for A = —0.3, and some gain 
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remains also for A = 0. 

With a 16 4 lattice a comparison is, as a matter of fact, impossible, simply because 
the separation between the phases is so strong that with conventional algorithms 
one does not observe any transition at all. With our algorithm, instead, on a 16 4 
lattice and for A = 0, we observe average tunneling times of the order of 10 3 (for 
tunneling time we follow the definition of fl3|). If we were interested only in A = 0, 
this number ought to be multiplied by 25, i.e. the number of X q involved. This 
is certainly not a small time, however it is small as compared to infinity, which 
corresponds to observing no transition at all. 

For a further reduction of the autocorrelation times, in addition to circumventing 
tunneling, one would have to replace the local Metropolis steps for with more 
efficient ones. In the dynamical parameter algorithm for the Potts model the 
cluster steps, which were originally introduced to implement the transitions between 
different q, have the additional advantage of reducing critical slowing down and, 
correspondingly, the autocorrelation time in the second order region. However, at 
this stage an implementation of cluster steps for gauge theories with continous groups 
appears very problematic, if not plainly impossible HT4| . A more promising direction 
to pursue might be along the lines of multi-scale algorithms | 15H , provided that these 
could be modified to account for the actual structure of the configurations. 
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Table 1 



(3(X) and g(X) of simulations on 16 4 lattice 
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Figure 1: Location of phase transition as function of A for 8 4 (circles) and 16 4 
(crosses) lattices. 




Figure 2: Distribution P(E, A) at the phase transition line on 16 4 lattice. 
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